Exploring a ZDR Offset Value
Contents
Exploring a ZDR Offset ValueΒΆ
import xarray as xr
import pyart
import numpy as np
import hvplot.xarray
import hvplot.pandas
import holoviews as hv
import warnings
import xcollection as xc
import glob
import matplotlib.pyplot as plt
import pandas as pd
import cartopy.crs as ccrs
from distributed import Client, LocalCluster
warnings.filterwarnings("ignore")
hv.extension('bokeh')
Read in the dataΒΆ
def create_dataset_from_sweep(radar, sweep=0, field=None):
"""
Creates an xarray.Dataset from sweeps
"""
# Grab the range
range = radar.range['data']
# Grab the elevation
elevation = radar.get_elevation(sweep)
# Grab the azimuth
azimuth = radar.get_azimuth(sweep)
# Grab the x, y, z values
x, y, z = radar.get_gate_x_y_z(sweep)
# Grab the lat, lon, and elevation
lat, lon, alt = radar.get_gate_lat_lon_alt(sweep)
# Add the fields
field_dict = {}
for field in list(radar.fields):
field_dict[field]=(["azimuth", "range"], radar.get_field(sweep, field))
ds = xr.Dataset(
data_vars=field_dict,
coords=dict(
x=(["azimuth", "range"], x),
y=(["azimuth", "range"], y),
z=(["azimuth", "range"], z),
lat=(["azimuth", "range"], lat),
lon=(["azimuth", "range"], lon),
alt=(["azimuth", "range"], alt),
azimuth=azimuth,
elevation=elevation,
sweep=np.array([sweep]),
range=range),
)
return ds.chunk()
def convert_to_xradar(radar):
"""
Converts from radar to xradar
"""
ds_list = []
sweeps = radar.sweep_number['data']
for sweep in sweeps:
ds_list.append(create_dataset_from_sweep(radar, sweep))
# Convert the numpy array to a list of strings
dict_keys = [x for x in list(sweeps.astype(str))]
# Zip the keys and datasets and turn into a dictionary
dict_of_dsets = dict(zip(dict_keys, ds_list))
return xc.Collection(dict_of_dsets)
List files availableΒΆ
We grab the files from the following directory, after ordering from the data portal and unzipping the archive
radar_files = sorted(glob.glob('../../data/x-band/ppi/*6_PPI.nc'))
Spin up a Dask ClusterΒΆ
cluster = LocalCluster()
client = Client(cluster)
client
distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-p297zaac', purging
distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-w0kfqa11', purging
distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-t_qx4ojw', purging
distributed.diskutils - INFO - Found stale lock file and directory '/Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-wpkibo54', purging
Client
Client-3d51261e-abc0-11ec-b6c4-acde48001122
| Connection method: Cluster object | Cluster type: distributed.LocalCluster |
| Dashboard: http://127.0.0.1:51851/status |
Cluster Info
LocalCluster
46ae0ffa
| Dashboard: http://127.0.0.1:51851/status | Workers: 4 |
| Total threads: 12 | Total memory: 16.00 GiB |
| Status: running | Using processes: True |
Scheduler Info
Scheduler
Scheduler-02e53ed5-9d8e-44b8-8b33-088e2998c32e
| Comm: tcp://127.0.0.1:51852 | Workers: 4 |
| Dashboard: http://127.0.0.1:51851/status | Total threads: 12 |
| Started: Just now | Total memory: 16.00 GiB |
Workers
Worker: 0
| Comm: tcp://127.0.0.1:51874 | Total threads: 3 |
| Dashboard: http://127.0.0.1:51875/status | Memory: 4.00 GiB |
| Nanny: tcp://127.0.0.1:51858 | |
| Local directory: /Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-x7jmkvn3 | |
Worker: 1
| Comm: tcp://127.0.0.1:51871 | Total threads: 3 |
| Dashboard: http://127.0.0.1:51872/status | Memory: 4.00 GiB |
| Nanny: tcp://127.0.0.1:51856 | |
| Local directory: /Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-x1bq8kb9 | |
Worker: 2
| Comm: tcp://127.0.0.1:51865 | Total threads: 3 |
| Dashboard: http://127.0.0.1:51867/status | Memory: 4.00 GiB |
| Nanny: tcp://127.0.0.1:51855 | |
| Local directory: /Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-k24izda2 | |
Worker: 3
| Comm: tcp://127.0.0.1:51866 | Total threads: 3 |
| Dashboard: http://127.0.0.1:51869/status | Memory: 4.00 GiB |
| Nanny: tcp://127.0.0.1:51857 | |
| Local directory: /Users/mgrover/git_repos/SAILS-radar-analysis/notebooks/march-14-2022/dask-worker-space/worker-dwu56dul | |
Loop through and read in the dataΒΆ
ds = convert_to_xradar(pyart.io.read(radar_files[0]))['0']
Clean the DataΒΆ
We will drop a couple of coordinates, and make sure that we donβt include the fill values (which are large negative numbers)
ds = ds.drop(['elevation', 'sweep'])
Set which fields to plot in the interactive plot
fields = ['DBZ', 'VEL', 'ZDR']
for field in fields:
da = ds[field]
ds[field] = ds[field].where(da > -999)
ds.compute().hvplot.quadmesh(x='lon',
y='lat',
ylabel='Latitude',
xlabel='Longitude',
z=['DBZ', 'VEL', 'ZDR'],
rasterize=True,
width=600,
height=400,
geo=True,
alpha=0.1,
selection_color={'NaN': (0, 0, 0, 0)},
#projection=ccrs.PlateCarree(),
cmap='pyart_Carbone42',
features={'rivers':'110m'},
tiles='StamenTerrainRetina').opts("Image", alpha=0.8)
Investigate HistogramsΒΆ
ReflectivityΒΆ
dbz_hist = ds.DBZ.hvplot.hist(bins=np.arange(-40, 40, 1), title='DBZ Values')
dbz_hist
NCPΒΆ
ncp_hist = ds.NCP.hvplot.hist(bins=np.arange(0, 1.01, 0.01), title='NCP Values')
ncp_hist
ZDRΒΆ
zdr_hist = ds.ZDR.hvplot.hist(bins=np.arange(-3, 1.1, .1), title='ZDR Values')
zdr_hist
Plot a 2D histogram using matplotlibΒΆ
ds_subset = ds.where((ds.RHOHV > 0.95) & (ds.NCP > 0.5))
bins = np.linspace(-15, 25, 71)
zdr_bins = np.linspace(-1, 0, 71)
hist = np.histogram2d(ds_subset.ZDR.values.flatten(), ds_subset.DBZ.values.flatten(), bins=70, range=((-1, 0), (-15, 25)))[0]
hist = np.ma.masked_where(hist == 0, hist)
# Create a 1-1 comparison
#x, y = np.meshgrid((-5, 1),(-6, 20))
c = plt.pcolormesh(bins,
zdr_bins,
hist,
cmap='pyart_HomeyerRainbow')
# Add a colorbar and labels
plt.colorbar(c, label='$log_{10}$ counts')
plt.xlabel('DBZ')
plt.ylabel('ZDR')
plt.show()
Plot a 2D Hexbin VisualizationΒΆ
zdr = ds_subset.ZDR.values.flatten()
dbz = ds_subset.DBZ.values.flatten()
df = pd.DataFrame({'ZDR':zdr,
'DBZ':dbz})
min_zdr = -2
max_zdr = 1
df_range = df[(df.ZDR >= min_zdr) &
(df.ZDR < max_zdr)]
df_range.hvplot.hexbin('DBZ',
'ZDR',
rasterize=True,
cmap='pyart_HomeyerRainbow',
logz=True,
)
Plot the Filtered DatasetΒΆ
def plot_field(field='DBZ',
color_bounds=(-20, 40),
cmap='pyart_Carbone42',
ds_to_plot=ds_subset):
plot = ds_to_plot[field].compute().hvplot.quadmesh(x='lon',
y='lat',
ylabel='Latitude',
xlabel='Longitude',
clim=color_bounds,
rasterize=True,
xlim=(-107.1, -106.8),
ylim=(38.8, 39.1),
width=600,
height=400,
geo=True,
alpha=0.1,
projection=ccrs.PlateCarree(),
cmap=cmap,
features={'rivers':'110m'},
tiles='StamenTerrainRetina').opts("Image",
alpha=0.8)
return plot
ref_plot = plot_field('DBZ')
vel_plot = plot_field('VEL', (-15, 15), cmap='pyart_balance')
zdr_plot = plot_field('ZDR', (-2, 1))
(ref_plot + vel_plot + zdr_plot).cols(1)